✅ Key points

  1. Opening line: starts with three backticks + {r ...}
    • Contains the chunk name (chunk_name) and options (echo, eval, etc.)
  2. Code: written below the opening line
  3. Closing line: three backticks alone (\```)
    • Nothing else should be on this line
  4. Optional: you can add comments inside the chunk (#) to describe what the code does
Option Purpose
echo=TRUE Shows the R code
eval=TRUE Runs the code (otherwise nothing is executed)
message=FALSE Hides package startup messages
warning=FALSE Hides warnings in the HTML
complete_2023_2024 <- read.csv("input/2023_2024_all_moth_counts.csv")
#clean_complete - full summer complete moth counts available for 2023 and 2024

library(tidyverse)
library(lme4)
library(performance)
library(summarytools)
library(ggeffects)
library(marginaleffects)
library(brms)
library(car)
library(viridis)
library(broom)
library(knitr)
library(dplyr)
library(kableExtra)
library(gt)
library(modelsummary)
library(sjPlot)

# Data Cleaning -----------------------------------------------------------


##To explore the distribution of your variables and count data like clean_complete
# quick visualizations
dfSummary(complete_2023_2024)
## Data Frame Summary  
## complete_2023_2024  
## Dimensions: 255 x 20  
## Duplicates: 0  
## 
## ------------------------------------------------------------------------------------------------------------------------
## No   Variable            Stats / Values                Freqs (% of Valid)    Graph                 Valid      Missing   
## ---- ------------------- ----------------------------- --------------------- --------------------- ---------- ----------
## 1    patch_name          1. Oka                        44 (17.3%)            III                   255        0         
##      [character]         2. Orford                     39 (15.3%)            III                   (100.0%)   (0.0%)    
##                          3. Rigaud                     30 (11.8%)            II                                         
##                          4. Mont_Saint_Hilaire         29 (11.4%)            II                                         
##                          5. Montebello                 22 ( 8.6%)            I                                          
##                          6. Notre_Dame_de_Bonsecours   18 ( 7.1%)            I                                          
##                          7. Mont_Royal                 13 ( 5.1%)            I                                          
##                          8. Mont_Saint_Bruno           12 ( 4.7%)                                                       
##                          9. Papineauville              12 ( 4.7%)                                                       
##                          10. Brownsburg                 6 ( 2.4%)                                                       
##                          [ 5 others ]                  30 (11.8%)            II                                         
## 
## 2    Years               Min  : 2023                   2023 : 102 (40.0%)    IIIIIIII              255        0         
##      [integer]           Mean : 2023.6                 2024 : 153 (60.0%)    IIIIIIIIIIII          (100.0%)   (0.0%)    
##                          Max  : 2024                                                                                    
## 
## 3    Year                Min  : 0                      0 : 102 (40.0%)       IIIIIIII              255        0         
##      [integer]           Mean : 0.6                    1 : 153 (60.0%)       IIIIIIIIIIII          (100.0%)   (0.0%)    
##                          Max  : 1                                                                                       
## 
## 4    stand_type          1. MOM                        10 ( 3.9%)                                  255        0         
##      [character]         2. Oak                        73 (28.6%)            IIIII                 (100.0%)   (0.0%)    
##                          3. Oak/Other                  33 (12.9%)            II                                         
##                          4. Oak/Pine                   41 (16.1%)            III                                        
##                          5. Other                      30 (11.8%)            II                                         
##                          6. Pine                       39 (15.3%)            III                                        
##                          7. Pine/Oak                   29 (11.4%)            II                                         
## 
## 5    landscape_type      1. agricultural               139 (54.5%)           IIIIIIIIII            255        0         
##      [character]         2. forest·                     85 (33.3%)           IIIIII                (100.0%)   (0.0%)    
##                          3. urban                       31 (12.2%)           II                                         
## 
## 6    stand_category      1. C                          35 (17.9%)            III                   196        59        
##      [character]         2. I                          29 (14.8%)            II                    (76.9%)    (23.1%)   
##                          3. E                          24 (12.2%)            II                                         
##                          4. H                          18 ( 9.2%)            I                                          
##                          5. F                          17 ( 8.7%)            I                                          
##                          6. L                          15 ( 7.7%)            I                                          
##                          7. A                          10 ( 5.1%)            I                                          
##                          8. D                          10 ( 5.1%)            I                                          
##                          9. K                          10 ( 5.1%)            I                                          
##                          10. MOM                       10 ( 5.1%)            I                                          
##                          [ 4 others ]                  18 ( 9.2%)            I                                          
## 
## 7    Percent_Oak         Mean (sd) : 0.3 (0.2)         0.00 : 48 (18.8%)     III                   255        0         
##      [numeric]           min < med < max:              0.10 : 22 ( 8.6%)     I                     (100.0%)   (0.0%)    
##                          0 < 0.3 < 0.8                 0.20 : 10 ( 3.9%)                                                
##                          IQR (CV) : 0.4 (0.7)          0.30 : 58 (22.7%)     IIII                                       
##                                                        0.40 : 50 (19.6%)     III                                        
##                                                        0.50 : 18 ( 7.1%)     I                                          
##                                                        0.60 : 18 ( 7.1%)     I                                          
##                                                        0.70 : 19 ( 7.5%)     I                                          
##                                                        0.80 : 12 ( 4.7%)                                                
## 
## 8    Percent_Pine        Mean (sd) : 0.2 (0.2)         0.00 : 117 (45.9%)    IIIIIIIII             255        0         
##      [numeric]           min < med < max:              0.10 :  29 (11.4%)    II                    (100.0%)   (0.0%)    
##                          0 < 0.1 < 0.8                 0.20 :  27 (10.6%)    II                                         
##                          IQR (CV) : 0.4 (1.2)          0.30 :  12 ( 4.7%)                                               
##                                                        0.40 :  21 ( 8.2%)    I                                          
##                                                        0.50 :  19 ( 7.5%)    I                                          
##                                                        0.60 :   9 ( 3.5%)                                               
##                                                        0.70 :  15 ( 5.9%)    I                                          
##                                                        0.80 :   6 ( 2.4%)                                               
## 
## 9    Percent_Maple       Mean (sd) : 0.3 (0.2)         10 distinct values    II                    246        9         
##      [numeric]           min < med < max:                                    III                   (96.5%)    (3.5%)    
##                          0 < 0.2 < 0.9                                       IIII                                       
##                          IQR (CV) : 0.3 (0.6)                                IIII                                       
##                                                                              IIII                                       
## 
## 10   Forest_Type         1. oak                        65 (25.5%)            IIIII                 255        0         
##      [character]         2. oak_maple                  48 (18.8%)            III                   (100.0%)   (0.0%)    
##                          3. pine                       45 (17.6%)            III                                        
##                          4. oak_pine                   21 ( 8.2%)            I                                          
##                          5. maple_oak                  16 ( 6.3%)            I                                          
##                          6. pine_oak                   14 ( 5.5%)            I                                          
##                          7. maple                      12 ( 4.7%)                                                       
##                          8. maple_birch                 7 ( 2.7%)                                                       
##                          9. maple_hardwood?             5 ( 2.0%)                                                       
##                          10. pine_hemlock               5 ( 2.0%)                                                       
##                          [ 9 others ]                  17 ( 6.7%)            I                                          
## 
## 11   stand_area_ha       Mean (sd) : 8.4 (3.5)         67 distinct values      . :                 255        0         
##      [numeric]           min < med < max:                                    : : :                 (100.0%)   (0.0%)    
##                          3.5 < 8 < 25.9                                      : : :                                      
##                          IQR (CV) : 3.2 (0.4)                                : : : :                                    
##                                                                              : : : : :   .                              
## 
## 12   forest_area_km2     Mean (sd) : 210.2 (436.8)     13 distinct values    :                     255        0         
##      [numeric]           min < med < max:                                    :                     (100.0%)   (0.0%)    
##                          1.6 < 14 < 1282                                     :                                          
##                          IQR (CV) : 95.7 (2.1)                               :                                          
##                                                                              :           :                              
## 
## 13   trap_name           1. BC Co-Dom1                   1 ( 0.4%)                                 255        0         
##      [character]         2. BC Co-Dom2                   1 ( 0.4%)                                 (100.0%)   (0.0%)    
##                          3. BC Dom1                      1 ( 0.4%)                                                      
##                          4. BC Dom2                      1 ( 0.4%)                                                      
##                          5. BC Low 1                     1 ( 0.4%)                                                      
##                          6. BC Low 2                     1 ( 0.4%)                                                      
##                          7. Bolt Co-Dom1                 1 ( 0.4%)                                                      
##                          8. Bolt Co-Dom2                 1 ( 0.4%)                                                      
##                          9. Bolt Dom1                    1 ( 0.4%)                                                      
##                          10. Bolt Dom2                   1 ( 0.4%)                                                      
##                          [ 245 others ]                245 (96.1%)           IIIIIIIIIIIIIIIIIII                        
## 
## 14   longitude           Mean (sd) : -73.7 (1)         255 distinct values       :                 255        0         
##      [numeric]           min < med < max:                                        :                 (100.0%)   (0.0%)    
##                          -75 < -74 < -72                                       : :   :   :                              
##                          IQR (CV) : 1.2 (0)                                    : :   :   :                              
##                                                                              : : : : :   : .                            
## 
## 15   total_moth_count    Mean (sd) : 57.6 (51.4)       113 distinct values   :                     252        3         
##      [integer]           min < med < max:                                    :                     (98.8%)    (1.2%)    
##                          0 < 47.5 < 378                                      : .                                        
##                          IQR (CV) : 55.5 (0.9)                               : :                                        
##                                                                              : : :                                      
## 
## 16   complete            Mean (sd) : 72.8 (60.9)       92 distinct values    :                     142        113       
##      [integer]           min < med < max:                                    : :                   (55.7%)    (44.3%)   
##                          3 < 61 < 378                                        : :                                        
##                          IQR (CV) : 68.8 (0.8)                               : : :                                      
##                                                                              : : : .   .                                
## 
## 17   clean_complete      Mean (sd) : 62.3 (55.4)       140 distinct values   :                     200        55        
##      [numeric]           min < med < max:                                    : .                   (78.4%)    (21.6%)   
##                          0 < 49.8 < 378                                      : :                                        
##                          IQR (CV) : 59 (0.9)                                 : :                                        
##                                                                              : : :                                      
## 
## 18   censored            Min  : 0                      0 : 200 (78.4%)       IIIIIIIIIIIIIII       255        0         
##      [integer]           Mean : 0.2                    1 :  55 (21.6%)       IIII                  (100.0%)   (0.0%)    
##                          Max  : 1                                                                                       
## 
## 19   reason_incomplete   1. (Empty string)             84 (60.4%)            IIIIIIIIIIII          139        116       
##      [character]         2. Missing in field            2 ( 1.4%)                                  (54.5%)    (45.5%)   
##                          3. Muck                       44 (31.7%)            IIIIII                                     
##                          4. Trap damaged·               9 ( 6.5%)            I                                          
## 
## 20   X                   All NA's                                                                  0          255       
##      [logical]                                                                                     (0.0%)     (100.0%)  
## ------------------------------------------------------------------------------------------------------------------------
str(complete_2023_2024)
## 'data.frame':    255 obs. of  20 variables:
##  $ patch_name       : chr  "Mont_Royal" "Mont_Royal" "Mont_Royal" "Mont_Royal" ...
##  $ Years            : int  2024 2024 2024 2024 2024 2024 2024 2024 2024 2024 ...
##  $ Year             : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ stand_type       : chr  "Oak" "Oak" "Oak" "Oak" ...
##  $ landscape_type   : chr  "urban" "urban" "urban" "urban" ...
##  $ stand_category   : chr  "E" "A" "C" "C" ...
##  $ Percent_Oak      : num  0.4 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.4 0.4 ...
##  $ Percent_Pine     : num  0 0 0 0 0 0 0 0 0.1 0.1 ...
##  $ Percent_Maple    : num  0.3 0.1 0.3 0.3 0.3 0.3 0.3 0.3 0 0 ...
##  $ Forest_Type      : chr  "oak_maple" "oak" "oak" "oak" ...
##  $ stand_area_ha    : num  14 6.4 8.6 8.6 8.6 11.7 11.7 11.7 8.9 8.9 ...
##  $ forest_area_km2  : num  4.89 4.89 4.89 4.89 4.89 4.89 4.89 4.89 4.89 4.89 ...
##  $ trap_name        : chr  "MRMOM" "MROak2.1" "MROak1.1" "MROak1.2" ...
##  $ longitude        : num  -73.6 -73.6 -73.6 -73.6 -73.6 ...
##  $ total_moth_count : int  93 15 7 20 6 96 52 NA 60 20 ...
##  $ complete         : int  93 15 7 20 6 96 52 NA 60 20 ...
##  $ clean_complete   : num  93 15 7 20 6 96 52 NA 60 20 ...
##  $ censored         : int  0 0 0 0 0 0 0 1 0 0 ...
##  $ reason_incomplete: chr  NA NA NA NA ...
##  $ X                : logi  NA NA NA NA NA NA ...
## change 'Year' from int to chr
complete_2023_2024$Year <- as.character(complete_2023_2024$Year)

str(complete_2023_2024)
## 'data.frame':    255 obs. of  20 variables:
##  $ patch_name       : chr  "Mont_Royal" "Mont_Royal" "Mont_Royal" "Mont_Royal" ...
##  $ Years            : int  2024 2024 2024 2024 2024 2024 2024 2024 2024 2024 ...
##  $ Year             : chr  "1" "1" "1" "1" ...
##  $ stand_type       : chr  "Oak" "Oak" "Oak" "Oak" ...
##  $ landscape_type   : chr  "urban" "urban" "urban" "urban" ...
##  $ stand_category   : chr  "E" "A" "C" "C" ...
##  $ Percent_Oak      : num  0.4 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.4 0.4 ...
##  $ Percent_Pine     : num  0 0 0 0 0 0 0 0 0.1 0.1 ...
##  $ Percent_Maple    : num  0.3 0.1 0.3 0.3 0.3 0.3 0.3 0.3 0 0 ...
##  $ Forest_Type      : chr  "oak_maple" "oak" "oak" "oak" ...
##  $ stand_area_ha    : num  14 6.4 8.6 8.6 8.6 11.7 11.7 11.7 8.9 8.9 ...
##  $ forest_area_km2  : num  4.89 4.89 4.89 4.89 4.89 4.89 4.89 4.89 4.89 4.89 ...
##  $ trap_name        : chr  "MRMOM" "MROak2.1" "MROak1.1" "MROak1.2" ...
##  $ longitude        : num  -73.6 -73.6 -73.6 -73.6 -73.6 ...
##  $ total_moth_count : int  93 15 7 20 6 96 52 NA 60 20 ...
##  $ complete         : int  93 15 7 20 6 96 52 NA 60 20 ...
##  $ clean_complete   : num  93 15 7 20 6 96 52 NA 60 20 ...
##  $ censored         : int  0 0 0 0 0 0 0 1 0 0 ...
##  $ reason_incomplete: chr  NA NA NA NA ...
##  $ X                : logi  NA NA NA NA NA NA ...
# remove space in trap_name
# replace '-' with '_' in trap_name
complete_2023_2024 <- complete_2023_2024 %>%
  mutate(trap_name = str_replace(trap_name, " ", ""))
complete_2023_2024 <- complete_2023_2024 %>%
  mutate(trap_name = str_replace(trap_name, "-", "_"))

#Remove MOM traps from either "stand type" or "stand category"
stand_type_filtered <- complete_2023_2024 %>%
  filter(stand_type != "MOM")
stand_category_filtered <- complete_2023_2024 %>%
  filter(stand_category != "MOM")



# looking for mistakes
unique(complete_2023_2024$stand_type)
## [1] "Oak"       "Oak/Pine"  "Pine"      "MOM"       "Pine/Oak"  "Oak/Other"
## [7] "Other"
unique(complete_2023_2024$patch_name)
##  [1] "Mont_Royal"               "Mont_Saint_Hilaire"      
##  [3] "Montebello"               "Notre_Dame_de_Bonsecours"
##  [5] "Oka"                      "Orford"                  
##  [7] "Papineauville"            "Rigaud"                  
##  [9] "Brownsburg"               "Hatley"                  
## [11] "Kenauk"                   "Mont_Gauvin/Glen "       
## [13] "Mont_Saint_Bruno"         "Parc_Michel_Chartrand"   
## [15] "Parc_Pointe_aux_Prairies"
unique(complete_2023_2024$Year)
## [1] "1" "0"
unique(complete_2023_2024$landscape_type)
## [1] "urban"        "agricultural" "forest "
unique(complete_2023_2024$stand_category)
##  [1] "E"   "A"   "C"   "L"   "MOM" "D"   "H"   "I"   "B"   "G"   "K"   "M"  
## [13] "F"   "J"   NA
unique(stand_type_filtered$stand_type)
## [1] "Oak"       "Oak/Pine"  "Pine"      "Pine/Oak"  "Oak/Other" "Other"
unique(stand_category_filtered$stand_category)
##  [1] "E" "A" "C" "L" "D" "H" "I" "B" "G" "K" "M" "F" "J"
# Data Summaries ----------------------------------------------------------

##separate Stand Type column so that we have a Stand ID for each stand in each patch
stand_ID_filtered <- stand_type_filtered %>% 
  separate(trap_name, into = c("stand_ID", "trap_number"), remove = FALSE, sep = "\\." ) %>% 
  glimpse()
## Rows: 245
## Columns: 22
## $ patch_name        <chr> "Mont_Royal", "Mont_Royal", "Mont_Royal", "Mont_Roya…
## $ Years             <int> 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2024…
## $ Year              <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1…
## $ stand_type        <chr> "Oak", "Oak", "Oak", "Oak", "Oak", "Oak", "Oak", "Oa…
## $ landscape_type    <chr> "urban", "urban", "urban", "urban", "urban", "urban"…
## $ stand_category    <chr> "E", "A", "C", "C", "C", "E", "E", "E", "E", "E", "E…
## $ Percent_Oak       <dbl> 0.4, 0.8, 0.6, 0.6, 0.6, 0.4, 0.4, 0.4, 0.4, 0.4, 0.…
## $ Percent_Pine      <dbl> 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1, 0.1, 0.…
## $ Percent_Maple     <dbl> 0.3, 0.1, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.0, 0.0, 0.…
## $ Forest_Type       <chr> "oak_maple", "oak", "oak", "oak", "oak", "oak_maple"…
## $ stand_area_ha     <dbl> 14.0, 6.4, 8.6, 8.6, 8.6, 11.7, 11.7, 11.7, 8.9, 8.9…
## $ forest_area_km2   <dbl> 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89…
## $ trap_name         <chr> "MRMOM", "MROak2.1", "MROak1.1", "MROak1.2", "MROak1…
## $ stand_ID          <chr> "MRMOM", "MROak2", "MROak1", "MROak1", "MROak1", "MR…
## $ trap_number       <chr> NA, "1", "1", "2", "3", "1", "2", "3", "1", "2", "3"…
## $ longitude         <dbl> -73.59218, -73.58859, -73.58863, -73.58738, -73.5892…
## $ total_moth_count  <int> 93, 15, 7, 20, 6, 96, 52, NA, 60, 20, 15, 17, 25, 3,…
## $ complete          <int> 93, 15, 7, 20, 6, 96, 52, NA, 60, 20, 15, 17, 25, 3,…
## $ clean_complete    <dbl> 93, 15, 7, 20, 6, 96, 52, NA, 60, 20, 15, 17, 25, 3,…
## $ censored          <int> 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1…
## $ reason_incomplete <chr> NA, NA, NA, NA, NA, NA, NA, "Missing in field", NA, …
## $ X                 <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
stand_ID_filtered_1 <- stand_category_filtered %>% 
  separate(trap_name, into = c("stand_ID", "trap_number"), remove = FALSE, sep = "\\." ) %>% 
  glimpse()
## Rows: 186
## Columns: 22
## $ patch_name        <chr> "Mont_Royal", "Mont_Royal", "Mont_Royal", "Mont_Roya…
## $ Years             <int> 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2024…
## $ Year              <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1…
## $ stand_type        <chr> "Oak", "Oak", "Oak", "Oak", "Oak", "Oak", "Oak", "Oa…
## $ landscape_type    <chr> "urban", "urban", "urban", "urban", "urban", "urban"…
## $ stand_category    <chr> "E", "A", "C", "C", "C", "E", "E", "E", "E", "E", "E…
## $ Percent_Oak       <dbl> 0.4, 0.8, 0.6, 0.6, 0.6, 0.4, 0.4, 0.4, 0.4, 0.4, 0.…
## $ Percent_Pine      <dbl> 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1, 0.1, 0.…
## $ Percent_Maple     <dbl> 0.3, 0.1, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.0, 0.0, 0.…
## $ Forest_Type       <chr> "oak_maple", "oak", "oak", "oak", "oak", "oak_maple"…
## $ stand_area_ha     <dbl> 14.0, 6.4, 8.6, 8.6, 8.6, 11.7, 11.7, 11.7, 8.9, 8.9…
## $ forest_area_km2   <dbl> 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89…
## $ trap_name         <chr> "MRMOM", "MROak2.1", "MROak1.1", "MROak1.2", "MROak1…
## $ stand_ID          <chr> "MRMOM", "MROak2", "MROak1", "MROak1", "MROak1", "MR…
## $ trap_number       <chr> NA, "1", "1", "2", "3", "1", "2", "3", "1", "2", "3"…
## $ longitude         <dbl> -73.59218, -73.58859, -73.58863, -73.58738, -73.5892…
## $ total_moth_count  <int> 93, 15, 7, 20, 6, 96, 52, NA, 60, 20, 15, 17, 25, 3,…
## $ complete          <int> 93, 15, 7, 20, 6, 96, 52, NA, 60, 20, 15, 17, 25, 3,…
## $ clean_complete    <dbl> 93, 15, 7, 20, 6, 96, 52, NA, 60, 20, 15, 17, 25, 3,…
## $ censored          <int> 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1…
## $ reason_incomplete <chr> NA, NA, NA, NA, NA, NA, NA, "Missing in field", NA, …
## $ X                 <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
total_moths_2023_2024 <- sum(complete_2023_2024$clean_complete, na.rm = TRUE)
print(total_moths_2023_2024)
## [1] 12456
n_distinct(complete_2023_2024$trap_name)
## [1] 255
n_distinct(stand_ID_filtered$stand_ID)
## [1] 155
n_distinct(stand_ID_filtered$patch_name)
## [1] 15
#check to see the distribution of moth count data
hist(complete_2023_2024$clean_complete, 
          main = " ", 
          xlab = "Spongy moth count/trap", 
          ylab = "Frequency", 
          col = "darkblue", 
          border = "black")

hist(complete_2023_2024$clean_complete, 
     main = " ", 
     xlab = "Spongy moth count/trap", 
     ylab = "Frequency", 
     col = "blue", 
     border = "black",
     breaks = 50)  # increase this number to make bins smaller

#Calculate the mean and standard deviation of moth counts for each 
#level of stand_category to see if there are differences.

# Checking unique combinations

table(stand_category_filtered$stand_category, 
      stand_category_filtered$patch_name)
##    
##     Brownsburg Hatley Kenauk Mont_Gauvin/Glen  Mont_Royal Mont_Saint_Bruno
##   A          0      0      1                 0          1                1
##   B          0      0      0                 0          0                0
##   C          2      2      1                 2          3                3
##   D          0      0      0                 0          0                0
##   E          0      0      0                 0          8                0
##   F          0      0      0                 0          0                0
##   G          0      0      0                 0          0                0
##   H          0      0      0                 0          0                0
##   I          0      0      0                 0          0                0
##   J          0      0      0                 0          0                0
##   K          0      0      0                 0          0                0
##   L          2      0      0                 0          1                0
##   M          0      0      0                 0          0                0
##    
##     Mont_Saint_Hilaire Montebello Notre_Dame_de_Bonsecours Oka Orford
##   A                  0          0                        0   0      3
##   B                  0          2                        0   0      0
##   C                  8          4                        6   0      2
##   D                  3          0                        0   0      0
##   E                  0          0                        0   8      6
##   F                  2          0                        5   0      5
##   G                  0          2                        0   1      0
##   H                  4          0                        1   8      0
##   I                  2          5                        0  18      0
##   J                  0          0                        0   0      4
##   K                  0          2                        0   5      0
##   L                  1          0                        0   2      7
##   M                  0          2                        2   0      2
##    
##     Papineauville Parc_Michel_Chartrand Rigaud
##   A             4                     0      0
##   B             0                     0      0
##   C             0                     2      0
##   D             0                     0      7
##   E             2                     0      0
##   F             0                     0      5
##   G             0                     0      0
##   H             2                     0      3
##   I             0                     0      4
##   J             0                     0      3
##   K             0                     0      3
##   L             0                     0      2
##   M             0                     0      0
table(stand_type_filtered$stand_type,
      stand_type_filtered$patch_name)
##            
##             Brownsburg Hatley Kenauk Mont_Gauvin/Glen  Mont_Royal
##   Oak                2      2      1                 2          8
##   Oak/Other          2      2      2                 2          0
##   Oak/Pine           0      0      1                 0          4
##   Other              0      2      2                 2          0
##   Pine               2      0      0                 0          1
##   Pine/Oak           0      0      0                 0          0
##            
##             Mont_Saint_Bruno Mont_Saint_Hilaire Montebello
##   Oak                      3                  8          6
##   Oak/Other                5                  3          2
##   Oak/Pine                 0                  8          2
##   Other                    4                  4          2
##   Pine                     0                  1          5
##   Pine/Oak                 0                  2          4
##            
##             Notre_Dame_de_Bonsecours Oka Orford Papineauville
##   Oak                              6   6     14             6
##   Oak/Other                        1   2      4             2
##   Oak/Pine                         5   9      2             2
##   Other                            2   0      4             2
##   Pine                             2   7     13             0
##   Pine/Oak                         1  18      0             0
##            
##             Parc_Michel_Chartrand Parc_Pointe_aux_Prairies Rigaud
##   Oak                           2                        0      7
##   Oak/Other                     2                        2      2
##   Oak/Pine                      0                        0      8
##   Other                         2                        4      0
##   Pine                          0                        0      8
##   Pine/Oak                      0                        0      4
# Set the levels of 'stand_type' to ensure the correct order
stand_type_filtered$stand_type <- factor(stand_type_filtered$stand_type, 
                              levels = c("Oak", "Oak/Pine", "Oak/Other", 
                                  "Pine/Oak", "Pine", "Other"))

# Create the 'stand-type by patch' table
table(stand_type_filtered$stand_type, stand_type_filtered$patch_name)
##            
##             Brownsburg Hatley Kenauk Mont_Gauvin/Glen  Mont_Royal
##   Oak                2      2      1                 2          8
##   Oak/Pine           0      0      1                 0          4
##   Oak/Other          2      2      2                 2          0
##   Pine/Oak           0      0      0                 0          0
##   Pine               2      0      0                 0          1
##   Other              0      2      2                 2          0
##            
##             Mont_Saint_Bruno Mont_Saint_Hilaire Montebello
##   Oak                      3                  8          6
##   Oak/Pine                 0                  8          2
##   Oak/Other                5                  3          2
##   Pine/Oak                 0                  2          4
##   Pine                     0                  1          5
##   Other                    4                  4          2
##            
##             Notre_Dame_de_Bonsecours Oka Orford Papineauville
##   Oak                              6   6     14             6
##   Oak/Pine                         5   9      2             2
##   Oak/Other                        1   2      4             2
##   Pine/Oak                         1  18      0             0
##   Pine                             2   7     13             0
##   Other                            2   0      4             2
##            
##             Parc_Michel_Chartrand Parc_Pointe_aux_Prairies Rigaud
##   Oak                           2                        0      7
##   Oak/Pine                      0                        0      8
##   Oak/Other                     2                        2      2
##   Pine/Oak                      0                        0      4
##   Pine                          0                        0      8
##   Other                         2                        4      0
# Create the contingency table
contingency_table <- table(stand_type_filtered$stand_type, stand_type_filtered$patch_name)

# Convert the table to a data frame
contingency_df <- as.data.frame(contingency_table)


# Heatmap stands by patch -------------------------------------------------

# Create a plot (heatmap) of the contingency table
contingency_df 
##         Var1                     Var2 Freq
## 1        Oak               Brownsburg    2
## 2   Oak/Pine               Brownsburg    0
## 3  Oak/Other               Brownsburg    2
## 4   Pine/Oak               Brownsburg    0
## 5       Pine               Brownsburg    2
## 6      Other               Brownsburg    0
## 7        Oak                   Hatley    2
## 8   Oak/Pine                   Hatley    0
## 9  Oak/Other                   Hatley    2
## 10  Pine/Oak                   Hatley    0
## 11      Pine                   Hatley    0
## 12     Other                   Hatley    2
## 13       Oak                   Kenauk    1
## 14  Oak/Pine                   Kenauk    1
## 15 Oak/Other                   Kenauk    2
## 16  Pine/Oak                   Kenauk    0
## 17      Pine                   Kenauk    0
## 18     Other                   Kenauk    2
## 19       Oak        Mont_Gauvin/Glen     2
## 20  Oak/Pine        Mont_Gauvin/Glen     0
## 21 Oak/Other        Mont_Gauvin/Glen     2
## 22  Pine/Oak        Mont_Gauvin/Glen     0
## 23      Pine        Mont_Gauvin/Glen     0
## 24     Other        Mont_Gauvin/Glen     2
## 25       Oak               Mont_Royal    8
## 26  Oak/Pine               Mont_Royal    4
## 27 Oak/Other               Mont_Royal    0
## 28  Pine/Oak               Mont_Royal    0
## 29      Pine               Mont_Royal    1
## 30     Other               Mont_Royal    0
## 31       Oak         Mont_Saint_Bruno    3
## 32  Oak/Pine         Mont_Saint_Bruno    0
## 33 Oak/Other         Mont_Saint_Bruno    5
## 34  Pine/Oak         Mont_Saint_Bruno    0
## 35      Pine         Mont_Saint_Bruno    0
## 36     Other         Mont_Saint_Bruno    4
## 37       Oak       Mont_Saint_Hilaire    8
## 38  Oak/Pine       Mont_Saint_Hilaire    8
## 39 Oak/Other       Mont_Saint_Hilaire    3
## 40  Pine/Oak       Mont_Saint_Hilaire    2
## 41      Pine       Mont_Saint_Hilaire    1
## 42     Other       Mont_Saint_Hilaire    4
## 43       Oak               Montebello    6
## 44  Oak/Pine               Montebello    2
## 45 Oak/Other               Montebello    2
## 46  Pine/Oak               Montebello    4
## 47      Pine               Montebello    5
## 48     Other               Montebello    2
## 49       Oak Notre_Dame_de_Bonsecours    6
## 50  Oak/Pine Notre_Dame_de_Bonsecours    5
## 51 Oak/Other Notre_Dame_de_Bonsecours    1
## 52  Pine/Oak Notre_Dame_de_Bonsecours    1
## 53      Pine Notre_Dame_de_Bonsecours    2
## 54     Other Notre_Dame_de_Bonsecours    2
## 55       Oak                      Oka    6
## 56  Oak/Pine                      Oka    9
## 57 Oak/Other                      Oka    2
## 58  Pine/Oak                      Oka   18
## 59      Pine                      Oka    7
## 60     Other                      Oka    0
## 61       Oak                   Orford   14
## 62  Oak/Pine                   Orford    2
## 63 Oak/Other                   Orford    4
## 64  Pine/Oak                   Orford    0
## 65      Pine                   Orford   13
## 66     Other                   Orford    4
## 67       Oak            Papineauville    6
## 68  Oak/Pine            Papineauville    2
## 69 Oak/Other            Papineauville    2
## 70  Pine/Oak            Papineauville    0
## 71      Pine            Papineauville    0
## 72     Other            Papineauville    2
## 73       Oak    Parc_Michel_Chartrand    2
## 74  Oak/Pine    Parc_Michel_Chartrand    0
## 75 Oak/Other    Parc_Michel_Chartrand    2
## 76  Pine/Oak    Parc_Michel_Chartrand    0
## 77      Pine    Parc_Michel_Chartrand    0
## 78     Other    Parc_Michel_Chartrand    2
## 79       Oak Parc_Pointe_aux_Prairies    0
## 80  Oak/Pine Parc_Pointe_aux_Prairies    0
## 81 Oak/Other Parc_Pointe_aux_Prairies    2
## 82  Pine/Oak Parc_Pointe_aux_Prairies    0
## 83      Pine Parc_Pointe_aux_Prairies    0
## 84     Other Parc_Pointe_aux_Prairies    4
## 85       Oak                   Rigaud    7
## 86  Oak/Pine                   Rigaud    8
## 87 Oak/Other                   Rigaud    2
## 88  Pine/Oak                   Rigaud    4
## 89      Pine                   Rigaud    8
## 90     Other                   Rigaud    0
ggplot(contingency_df, aes(x = Var1, y = Var2, fill = Freq)) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "darkred") +
  labs(x = "Stand Type", y = "Patch Name", fill = "Frequency") +
  theme_minimal()

# Save the plot as an image (e.g., PNG)
#ggsave("contingency_table_heatmap.png", width = 7.5, height = 6)

#convert the data in the heatmap plot to a table format
wide_table <- contingency_df %>%
  pivot_wider(names_from = Var1, values_from = Freq, values_fill = 0)

wide_table <- wide_table %>%
  mutate(
    Var2 = trimws(as.character(Var2)),
    Var2 = dplyr::recode(Var2,
                          "Mont_Gauvin/Glen" = "Mont Gauvin",
                          "Mont_Royal" = "Mont Royal",
                          "Mont_Saint_Bruno" = "Mont Saint Bruno",
                          "Mont_Saint_Hilaire" = "Mont Saint Hilaire",
                          "Notre_Dame_de_Bonsecours" = "Notre Dame de Bonsecours",
                          "Oka" = "Parc Oka",
                          "Orford" = "Mont Orford",
                          "Parc_Michel_Chartrand" = "Parc Michel Chartrand",
                          "Parc_Pointe_aux_Prairies" = "Parc Pointe aux Prairies",
                          "Rigaud" = "Mont Rigaud",
                         ))


##print(wide_table)

patch_order <- c("Papineauville", "Montebello", "Notre Dame de Bonsecours",
                 "Kenauk", "Brownsburg", "Mont Rigaud", "Parc Oka",
                 "Mont Royal",
                 "Parc Pointe aux Prairies", "Parc Michel Chartrand",
                 "Mont Saint Bruno", "Mont Saint Hilaire","Mont Gauvin",
                 "Mont Orford", "Hatley")
wide_table <- wide_table %>%
  mutate(Var2 = factor(Var2, levels = patch_order)) %>%
  arrange(Var2)

##print(wide_table)

gt_table <- wide_table %>%
  gt() %>%
  cols_label(
    Var2 = "Patch Name"
  ) %>%
  tab_options(
    table.font.size = 12,
    heading.title.font.size = 16,
    heading.subtitle.font.size = 14
  ) %>%
  cols_align(
    align = "center",
    columns = everything()
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels(everything())
  ) %>%
  tab_style(  # Custom striping on even-numbered rows
    style = cell_fill(color = "#f9f9f0"),  # You can change this color
    locations = cells_body(rows = seq(2, nrow(wide_table), 2))
  )


version$version.string
## [1] "R version 4.4.3 (2025-02-28)"
# Summary statistics for moth count by stand_category and stand_type
summary_stats <- stand_type_filtered %>%
  group_by(stand_category, stand_type, patch_name) %>%
  summarise(
    mean_count = mean(total_moth_count, na.rm = TRUE),
    sd_count = sd(total_moth_count, na.rm = TRUE),
    var_count = var(total_moth_count, na.rm = TRUE),
    count = n()
  )

##print(summary_stats, n=22)





# Table to use -----------------------------------------------------------

# Summary statistics for moth count by stand_type, differentiating between
#traps set and traps with usable data
summary_stats_3 <- stand_type_filtered %>%
  group_by(stand_type) %>%
  summarise(
    mean_count = mean(clean_complete, na.rm = TRUE),
    sd_count = sd(clean_complete, na.rm = TRUE),
    n_traps = n(),
    n_obs = sum(!is.na (clean_complete))
  )

print(summary_stats_3, n=22)
## # A tibble: 6 × 5
##   stand_type mean_count sd_count n_traps n_obs
##   <fct>           <dbl>    <dbl>   <int> <int>
## 1 Oak              77.0     69.8      73    57
## 2 Oak/Pine         64.4     52.6      41    30
## 3 Oak/Other        37.3     27.0      33    29
## 4 Pine/Oak         85.8     61.1      29    20
## 5 Pine             54.8     33.2      39    27
## 6 Other            36.9     26.8      30    27
# Round numeric columns to 2 decimal places
summary_stats_3$mean_count <- round(summary_stats_3$mean_count, 2)
summary_stats_3$sd_count <- round(summary_stats_3$sd_count, 2)
summary_stats_3$count <- round(summary_stats_3$n_obs, 2)



# Summary statistics for moth count by patch
summary_stats_4 <- stand_type_filtered %>%
  group_by(patch_name) %>%
  summarise(
    mean_count = mean(clean_complete, na.rm = TRUE),
    sd_count = sd(clean_complete, na.rm = TRUE),
    count = n()
  )

print(summary_stats_4, n=22)
## # A tibble: 15 × 4
##    patch_name                 mean_count sd_count count
##    <chr>                           <dbl>    <dbl> <int>
##  1 "Brownsburg"                     38.6     9.40     6
##  2 "Hatley"                         74.8    18.4      6
##  3 "Kenauk"                         26.8    20.4      6
##  4 "Mont_Gauvin/Glen "              34.1    11.6      6
##  5 "Mont_Royal"                     35.5    32.0     13
##  6 "Mont_Saint_Bruno"               24.7    13.9     12
##  7 "Mont_Saint_Hilaire"             23.5    29.6     26
##  8 "Montebello"                     54      44.8     21
##  9 "Notre_Dame_de_Bonsecours"       54.8    32.9     17
## 10 "Oka"                           105.     63.6     42
## 11 "Orford"                         88.0    67.9     37
## 12 "Papineauville"                  71.0    56.2     12
## 13 "Parc_Michel_Chartrand"          70.2    17.4      6
## 14 "Parc_Pointe_aux_Prairies"       43.3    32.7      6
## 15 "Rigaud"                         38.1    26.4     29
# Summary statistics for moth count by stand_type and patch
# summary_stats_5 <- stand_type_filtered %>%

# Remove MOM row in 'complete' moth counts
moth_by_stand_summary_stats <- summary_stats_3[-c(1),]
print(moth_by_stand_summary_stats, n=22)
## # A tibble: 5 × 6
##   stand_type mean_count sd_count n_traps n_obs count
##   <fct>           <dbl>    <dbl>   <int> <int> <dbl>
## 1 Oak/Pine         64.4     52.6      41    30    30
## 2 Oak/Other        37.3     27.0      33    29    29
## 3 Pine/Oak         85.8     61.1      29    20    20
## 4 Pine             54.8     33.2      39    27    27
## 5 Other            36.9     26.8      30    27    27
# Remove MOM row in 'clean_complete' moth counts
moth_by_stand_summary_stats_2 <- summary_stats_4[-c(1),]
print(moth_by_stand_summary_stats_2, n=22)
## # A tibble: 14 × 4
##    patch_name                 mean_count sd_count count
##    <chr>                           <dbl>    <dbl> <int>
##  1 "Hatley"                         74.8     18.4     6
##  2 "Kenauk"                         26.8     20.4     6
##  3 "Mont_Gauvin/Glen "              34.1     11.6     6
##  4 "Mont_Royal"                     35.5     32.0    13
##  5 "Mont_Saint_Bruno"               24.7     13.9    12
##  6 "Mont_Saint_Hilaire"             23.5     29.6    26
##  7 "Montebello"                     54       44.8    21
##  8 "Notre_Dame_de_Bonsecours"       54.8     32.9    17
##  9 "Oka"                           105.      63.6    42
## 10 "Orford"                         88.0     67.9    37
## 11 "Papineauville"                  71.0     56.2    12
## 12 "Parc_Michel_Chartrand"          70.2     17.4     6
## 13 "Parc_Pointe_aux_Prairies"       43.3     32.7     6
## 14 "Rigaud"                         38.1     26.4    29
# Visualizations ----------------------------------------------------------

##separate Stand Type column so that we have a Stand ID for each stand in each patch
#stand_ID_filtered <- stand_type_filtered %>% 

#Poisson, using all levels of data collection as a random effect, except Oak
#which is being fitted as a fixed effect
model_complete_poisson_oak <- glmer(
  round(clean_complete) ~ (1|trap_name) +
    (Percent_Oak) + (1|patch_name),
  family =poisson(), data = stand_ID_filtered)
summary(model_complete_poisson_oak)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(clean_complete) ~ (1 | trap_name) + (Percent_Oak) + (1 |  
##     patch_name)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1912.6    1925.6    -952.3    1904.6       186 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.06666 -0.11714  0.00811  0.08998  0.36154 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  trap_name  (Intercept) 0.588    0.7668  
##  patch_name (Intercept) 0.262    0.5118  
## Number of obs: 190, groups:  trap_name, 190; patch_name, 15
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.3637     0.1752  19.201  < 2e-16 ***
## Percent_Oak   0.7238     0.2630   2.752  0.00592 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Percent_Oak -0.525
performance::check_overdispersion(model_complete_poisson_oak)
## # Overdispersion test
## 
##        dispersion ratio =  0.142
##   Pearson's Chi-Squared = 26.378
##                 p-value =      1
performance::check_model(model_complete_poisson_oak)

#Heavily underdispersed, indicating that the moth counts are less variable
#than expected

#Run the same basic model, but have stands nested within patches as a 
#random effects, first keeping in trap_name as a random effect and then 
#removing it

model_complete_poisson_oak_nested <- glmer(
  round(clean_complete) ~ Percent_Oak + 
    (1 | trap_name) + 
    (1 | patch_name/stand_ID),   # nesting structure
  family = poisson(),
  data = stand_ID_filtered
)
summary(model_complete_poisson_oak_nested)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(clean_complete) ~ Percent_Oak + (1 | trap_name) + (1 |  
##     patch_name/stand_ID)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1905.2    1921.5    -947.6    1895.2       185 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.96329 -0.11714  0.00997  0.09144  0.37052 
## 
## Random effects:
##  Groups              Name        Variance Std.Dev.
##  trap_name           (Intercept) 0.3346   0.5785  
##  stand_ID:patch_name (Intercept) 0.2975   0.5454  
##  patch_name          (Intercept) 0.2328   0.4824  
## Number of obs: 190, groups:  
## trap_name, 190; stand_ID:patch_name, 131; patch_name, 15
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.3662     0.1785  18.862   <2e-16 ***
## Percent_Oak   0.6520     0.2997   2.175   0.0296 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Percent_Oak -0.578
performance::check_overdispersion(model_complete_poisson_oak_nested)
## # Overdispersion test
## 
##        dispersion ratio =  0.146
##   Pearson's Chi-Squared = 27.003
##                 p-value =      1
performance::check_model(model_complete_poisson_oak_nested)

model_complete_poisson_oak_nested_1 <- glmer(
  round(clean_complete) ~ Percent_Oak + 
    #(1 | trap_name) + 
    (1 | patch_name/stand_ID),   # nesting structure
  family = poisson(),
  data = stand_ID_filtered
)
summary(model_complete_poisson_oak_nested_1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(clean_complete) ~ Percent_Oak + (1 | patch_name/stand_ID)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2887.1    2900.1   -1439.6    2879.1       186 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.0329 -0.7169 -0.0228  0.1485 13.3173 
## 
## Random effects:
##  Groups              Name        Variance Std.Dev.
##  stand_ID:patch_name (Intercept) 0.6051   0.7779  
##  patch_name          (Intercept) 0.2183   0.4672  
## Number of obs: 190, groups:  stand_ID:patch_name, 131; patch_name, 15
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.3726     0.1777  18.982   <2e-16 ***
## Percent_Oak   0.6415     0.3083   2.081   0.0374 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Percent_Oak -0.596
performance::check_overdispersion(model_complete_poisson_oak_nested_1)
## # Overdispersion test
## 
##        dispersion ratio =    6.476
##   Pearson's Chi-Squared = 1204.591
##                 p-value =  < 0.001
performance::check_model(model_complete_poisson_oak_nested_1)

#Try 3 levels of random effects
model_complete_poisson_oak_nested_2 <- glmer(
  round(clean_complete) ~ Percent_Oak + 
    (1 | patch_name/stand_ID/trap_name),   # nesting structure
  family = poisson(),
  data = stand_ID_filtered
)
summary(model_complete_poisson_oak_nested_2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## round(clean_complete) ~ Percent_Oak + (1 | patch_name/stand_ID/trap_name)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1905.2    1921.5    -947.6    1895.2       185 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.96329 -0.11714  0.00997  0.09144  0.37052 
## 
## Random effects:
##  Groups                        Name        Variance Std.Dev.
##  trap_name:stand_ID:patch_name (Intercept) 0.3346   0.5785  
##  stand_ID:patch_name           (Intercept) 0.2975   0.5454  
##  patch_name                    (Intercept) 0.2328   0.4824  
## Number of obs: 190, groups:  
## trap_name:stand_ID:patch_name, 190; stand_ID:patch_name, 131; patch_name, 15
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.3662     0.1785  18.862   <2e-16 ***
## Percent_Oak   0.6520     0.2997   2.175   0.0296 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Percent_Oak -0.578
performance::check_overdispersion(model_complete_poisson_oak_nested_2)
## # Overdispersion test
## 
##        dispersion ratio =  0.146
##   Pearson's Chi-Squared = 27.003
##                 p-value =      1
performance::check_model(model_complete_poisson_oak_nested_2)

#Poisson, using all levels of data collection as a random effect, except Pine
#which is being fitted as a fixed effect
model_complete_poisson_pine <- glmer(
  round(clean_complete) ~ (1|trap_name) +
    (Percent_Pine) + (1|patch_name),
  family =poisson(), data = stand_ID_filtered)
summary(model_complete_poisson_pine)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(clean_complete) ~ (1 | trap_name) + (Percent_Pine) + (1 |  
##     patch_name)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1919.8    1932.8    -955.9    1911.8       186 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.94420 -0.11518  0.00723  0.08474  0.29169 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  trap_name  (Intercept) 0.6211   0.7881  
##  patch_name (Intercept) 0.2417   0.4916  
## Number of obs: 190, groups:  trap_name, 190; patch_name, 15
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.5990     0.1499  24.015   <2e-16 ***
## Percent_Pine   0.1201     0.2866   0.419    0.675    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Percent_Pin -0.245
performance::check_overdispersion(model_complete_poisson_pine)
## # Overdispersion test
## 
##        dispersion ratio =  0.137
##   Pearson's Chi-Squared = 25.519
##                 p-value =      1
performance::check_model(model_complete_poisson_pine)

#Again, heavily underdispersed, indicating that the moth counts are 
#less variable than expected
##Run the same basic model, but have stands nested within patches as a 
#random effects, first keeping in trap_name as a random effect and then 
#removing it

model_complete_poisson_pine_nested <- glmer(
  round(clean_complete) ~ (Percent_Pine) + 
    (1 | trap_name) + 
    (1 | patch_name/stand_ID),   # nesting structure
  family = poisson(),
  data = stand_ID_filtered)
summary(model_complete_poisson_pine_nested)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(clean_complete) ~ (Percent_Pine) + (1 | trap_name) + (1 |  
##     patch_name/stand_ID)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1908.4    1924.6    -949.2    1898.4       185 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.83251 -0.13372  0.00920  0.08245  0.39874 
## 
## Random effects:
##  Groups              Name        Variance Std.Dev.
##  trap_name           (Intercept) 0.3234   0.5687  
##  stand_ID:patch_name (Intercept) 0.3452   0.5875  
##  patch_name          (Intercept) 0.2095   0.4577  
## Number of obs: 190, groups:  
## trap_name, 190; stand_ID:patch_name, 131; patch_name, 15
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.5375     0.1478  23.936   <2e-16 ***
## Percent_Pine   0.4124     0.3502   1.177    0.239    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Percent_Pin -0.291
performance::check_overdispersion(model_complete_poisson_pine_nested)
## # Overdispersion test
## 
##        dispersion ratio =  0.142
##   Pearson's Chi-Squared = 26.303
##                 p-value =      1
performance::check_model(model_complete_poisson_pine_nested)

model_complete_poisson_pine_nested_1 <- glmer(
  round(clean_complete) ~ (Percent_Pine) + 
    #(1 | trap_name) + 
    (1 | patch_name/stand_ID),   # nesting structure
  family = poisson(),
  data = stand_ID_filtered)
summary(model_complete_poisson_pine_nested_1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(clean_complete) ~ (Percent_Pine) + (1 | patch_name/stand_ID)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2889.1    2902.1   -1440.6    2881.1       186 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.0307 -0.6664 -0.0092  0.1470 13.3232 
## 
## Random effects:
##  Groups              Name        Variance Std.Dev.
##  stand_ID:patch_name (Intercept) 0.6281   0.7925  
##  patch_name          (Intercept) 0.1994   0.4466  
## Number of obs: 190, groups:  stand_ID:patch_name, 131; patch_name, 15
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.5296     0.1451  24.331   <2e-16 ***
## Percent_Pine   0.5442     0.3650   1.491    0.136    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Percent_Pin -0.287
performance::check_overdispersion(model_complete_poisson_pine_nested_1)
## # Overdispersion test
## 
##        dispersion ratio =    6.474
##   Pearson's Chi-Squared = 1204.143
##                 p-value =  < 0.001
performance::check_model(model_complete_poisson_pine_nested_1)

model_complete_poisson_pine_nested_2 <- glmer(
  round(clean_complete) ~ (Percent_Pine) + 
    (1 | patch_name/stand_ID/trap_name),   # nesting structure
  family = poisson(),
  data = stand_ID_filtered)
summary(model_complete_poisson_pine_nested_2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## round(clean_complete) ~ (Percent_Pine) + (1 | patch_name/stand_ID/trap_name)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1908.4    1924.6    -949.2    1898.4       185 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.83251 -0.13372  0.00920  0.08245  0.39874 
## 
## Random effects:
##  Groups                        Name        Variance Std.Dev.
##  trap_name:stand_ID:patch_name (Intercept) 0.3234   0.5687  
##  stand_ID:patch_name           (Intercept) 0.3452   0.5875  
##  patch_name                    (Intercept) 0.2095   0.4577  
## Number of obs: 190, groups:  
## trap_name:stand_ID:patch_name, 190; stand_ID:patch_name, 131; patch_name, 15
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.5375     0.1478  23.936   <2e-16 ***
## Percent_Pine   0.4124     0.3502   1.177    0.239    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Percent_Pin -0.291
performance::check_overdispersion(model_complete_poisson_pine_nested_2)
## # Overdispersion test
## 
##        dispersion ratio =  0.142
##   Pearson's Chi-Squared = 26.303
##                 p-value =      1
performance::check_model(model_complete_poisson_pine_nested_2)

#model for Oak and Pine together
model_both <- glmer(round(clean_complete) ~ Percent_Pine + Percent_Oak + 
                      (1 | trap_name) + (1 | patch_name), 
                    data = stand_ID_filtered, family = poisson)

summary(model_both)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## round(clean_complete) ~ Percent_Pine + Percent_Oak + (1 | trap_name) +  
##     (1 | patch_name)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1909.4    1925.6    -949.7    1899.4       185 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.06551 -0.12913  0.00495  0.09740  0.32161 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  trap_name  (Intercept) 0.5761   0.7590  
##  patch_name (Intercept) 0.2482   0.4981  
## Number of obs: 190, groups:  trap_name, 190; patch_name, 15
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.1338     0.1995  15.705  < 2e-16 ***
## Percent_Pine   0.7514     0.3292   2.283 0.022448 *  
## Percent_Oak    1.1112     0.3113   3.569 0.000358 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Prcn_P
## Percent_Pin -0.509       
## Percent_Oak -0.660  0.546
performance::check_overdispersion(model_both)
## # Overdispersion test
## 
##        dispersion ratio =  0.141
##   Pearson's Chi-Squared = 26.063
##                 p-value =      1
#Again, heavily underdispersed, indicating that the moth counts are 
#less variable than expected
##Run the same basic model, but have stands nested within patches as a 
#random effects, first keeping in trap_name as a random effect and then 
#removing it

model_both_nested <- glmer(
            round(clean_complete) ~ Percent_Pine + Percent_Oak + 
            (1 | trap_name) + 
            (1 | patch_name/stand_ID), # nesting structure
            family = poisson(), 
            data = stand_ID_filtered)

summary(model_both_nested)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## round(clean_complete) ~ Percent_Pine + Percent_Oak + (1 | trap_name) +  
##     (1 | patch_name/stand_ID)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1901.0    1920.5    -944.5    1889.0       184 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.96859 -0.13088  0.00960  0.09818  0.38032 
## 
## Random effects:
##  Groups              Name        Variance Std.Dev.
##  trap_name           (Intercept) 0.3268   0.5716  
##  stand_ID:patch_name (Intercept) 0.2859   0.5347  
##  patch_name          (Intercept) 0.2239   0.4732  
## Number of obs: 190, groups:  
## trap_name, 190; stand_ID:patch_name, 131; patch_name, 15
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.1157     0.2023  15.399  < 2e-16 ***
## Percent_Pine   0.9440     0.3774   2.501  0.01238 *  
## Percent_Oak    1.0561     0.3368   3.136  0.00171 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Prcn_P
## Percent_Pin -0.498       
## Percent_Oak -0.678  0.473
performance::check_overdispersion(model_both_nested)
## # Overdispersion test
## 
##        dispersion ratio =  0.146
##   Pearson's Chi-Squared = 26.852
##                 p-value =      1
model_both_nested_1 <- glmer(
  round(clean_complete) ~ Percent_Pine + Percent_Oak + 
    #(1 | trap_name) + 
    (1 | patch_name/stand_ID), # nesting structure
  family = poisson(), 
  data = stand_ID_filtered)

summary(model_both_nested_1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## round(clean_complete) ~ Percent_Pine + Percent_Oak + (1 | patch_name/stand_ID)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    2881.9    2898.1   -1435.9    2871.9       185 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.0299 -0.7789 -0.0275  0.1466 13.3255 
## 
## Random effects:
##  Groups              Name        Variance Std.Dev.
##  stand_ID:patch_name (Intercept) 0.5731   0.7570  
##  patch_name          (Intercept) 0.2117   0.4601  
## Number of obs: 190, groups:  stand_ID:patch_name, 131; patch_name, 15
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.1210     0.1980  15.761  < 2e-16 ***
## Percent_Pine   1.0561     0.3896   2.711  0.00671 ** 
## Percent_Oak    1.0370     0.3349   3.097  0.00196 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Prcn_P
## Percent_Pin -0.476       
## Percent_Oak -0.677  0.438
performance::check_overdispersion(model_both_nested_1)
## # Overdispersion test
## 
##        dispersion ratio =    6.509
##   Pearson's Chi-Squared = 1204.198
##                 p-value =  < 0.001
model_both_nested_2 <- glmer(
  round(clean_complete) ~ Percent_Pine + Percent_Oak + 
    (1 | patch_name/stand_ID/trap_name), # nesting structure
  family = poisson(), 
  data = stand_ID_filtered)

summary(model_both_nested_2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## round(clean_complete) ~ Percent_Pine + Percent_Oak + (1 | patch_name/stand_ID/trap_name)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1901.0    1920.5    -944.5    1889.0       184 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.96859 -0.13088  0.00960  0.09818  0.38032 
## 
## Random effects:
##  Groups                        Name        Variance Std.Dev.
##  trap_name:stand_ID:patch_name (Intercept) 0.3268   0.5716  
##  stand_ID:patch_name           (Intercept) 0.2859   0.5347  
##  patch_name                    (Intercept) 0.2239   0.4732  
## Number of obs: 190, groups:  
## trap_name:stand_ID:patch_name, 190; stand_ID:patch_name, 131; patch_name, 15
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.1157     0.2023  15.399  < 2e-16 ***
## Percent_Pine   0.9440     0.3774   2.501  0.01238 *  
## Percent_Oak    1.0561     0.3368   3.136  0.00171 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Prcn_P
## Percent_Pin -0.498       
## Percent_Oak -0.678  0.473
performance::check_overdispersion(model_both_nested_2)
## # Overdispersion test
## 
##        dispersion ratio =  0.146
##   Pearson's Chi-Squared = 26.852
##                 p-value =      1
#model for Oak and Pine together, adding an interaction of oak & pine
model_both_2 <- glmer(round(clean_complete) ~ Percent_Pine + Percent_Oak + 
                      (1 | trap_name) + (1 | patch_name) + 
                      (Percent_Pine * Percent_Oak), 
                    data = stand_ID_filtered, family = poisson)

summary(model_both_2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## round(clean_complete) ~ Percent_Pine + Percent_Oak + (1 | trap_name) +  
##     (1 | patch_name) + (Percent_Pine * Percent_Oak)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1908.7    1928.2    -948.4    1896.7       184 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.04594 -0.13354  0.01572  0.09499  0.35363 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  trap_name  (Intercept) 0.5692   0.7544  
##  patch_name (Intercept) 0.2484   0.4984  
## Number of obs: 190, groups:  trap_name, 190; patch_name, 15
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                3.1521     0.1993  15.818  < 2e-16 ***
## Percent_Pine               0.4895     0.3640   1.345  0.17863    
## Percent_Oak                0.9508     0.3249   2.926  0.00343 ** 
## Percent_Pine:Percent_Oak   2.7525     1.6875   1.631  0.10286    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Prcn_P Prcn_O
## Percent_Pin -0.479              
## Percent_Oak -0.643  0.601       
## Prcnt_P:P_O  0.051 -0.437 -0.299
performance::check_overdispersion(model_both_2)
## # Overdispersion test
## 
##        dispersion ratio =  0.140
##   Pearson's Chi-Squared = 25.805
##                 p-value =      1
ggplot(stand_ID_filtered, aes(x = Percent_Oak, y = clean_complete)) +
  geom_point(alpha = 0.5) +
  stat_smooth(method = "glm", method.args = list(family = "poisson"), 
              se = TRUE, color = "darkgreen") +
  labs(
    title = "Effect of Percent Oak on Moth Counts",
    x = "% Oak Cover",
    y = "Moth Count"
  ) +
  theme_minimal()

# Checking Maple ----------------------------------------------------------

#Poisson, using all levels of data collection as a random effect, except Maple
#which is being fitted as a fixed effect
model_complete_poisson_maple <- glmer(
  round(clean_complete) ~ (1|trap_name) +
    (Percent_Maple) + (1|patch_name),
  family =poisson(), data = stand_ID_filtered)
summary(model_complete_poisson_maple)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(clean_complete) ~ (1 | trap_name) + (Percent_Maple) + (1 |  
##     patch_name)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1872.4    1885.3    -932.2    1864.4       182 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92817 -0.13025  0.00652  0.09274  0.30691 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  trap_name  (Intercept) 0.5817   0.7627  
##  patch_name (Intercept) 0.2590   0.5090  
## Number of obs: 186, groups:  trap_name, 186; patch_name, 15
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     3.7646     0.1869  20.144   <2e-16 ***
## Percent_Maple  -0.4860     0.4528  -1.073    0.283    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Percent_Mpl -0.604
performance::check_overdispersion(model_complete_poisson_maple)
## # Overdispersion test
## 
##        dispersion ratio =  0.127
##   Pearson's Chi-Squared = 23.172
##                 p-value =      1
performance::check_model(model_complete_poisson_maple)

#model for Oak and Maple together
model_oak_maple <- glmer(round(clean_complete) ~ Percent_Maple + Percent_Oak + 
                      (1 | trap_name) + (1 | patch_name), 
                    data = stand_ID_filtered, family = poisson)

summary(model_oak_maple)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## round(clean_complete) ~ Percent_Maple + Percent_Oak + (1 | trap_name) +  
##     (1 | patch_name)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1863.7    1879.8    -926.8    1853.7       181 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.05891 -0.12810  0.00800  0.08859  0.38383 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  trap_name  (Intercept) 0.5399   0.7348  
##  patch_name (Intercept) 0.2856   0.5344  
## Number of obs: 186, groups:  trap_name, 186; patch_name, 15
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     3.4304     0.2140  16.028  < 2e-16 ***
## Percent_Maple  -0.3332     0.4421  -0.754 0.451113    
## Percent_Oak     0.8734     0.2622   3.331 0.000864 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Prcn_M
## Percent_Mpl -0.556       
## Percent_Oak -0.470  0.098
performance::check_overdispersion(model_oak_maple)
## # Overdispersion test
## 
##        dispersion ratio =  0.132
##   Pearson's Chi-Squared = 23.921
##                 p-value =      1
#model for Oak and Maple together, adding an interaction of oak & maple
model_oak_maple_2 <- glmer(round(clean_complete) ~ Percent_Maple + 
                             Percent_Oak + 
                        (1 | trap_name) + (1 | patch_name) + 
                        (Percent_Maple * Percent_Oak), 
                      data = stand_ID_filtered, family = poisson)

summary(model_oak_maple_2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## round(clean_complete) ~ Percent_Maple + Percent_Oak + (1 | trap_name) +  
##     (1 | patch_name) + (Percent_Maple * Percent_Oak)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1865.6    1884.9    -926.8    1853.6       180 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.06075 -0.12872  0.00638  0.08843  0.37974 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  trap_name  (Intercept) 0.5401   0.7349  
##  patch_name (Intercept) 0.2842   0.5331  
## Number of obs: 186, groups:  trap_name, 186; patch_name, 15
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 3.4649     0.2400  14.440   <2e-16 ***
## Percent_Maple              -0.4672     0.6119  -0.763   0.4452    
## Percent_Oak                 0.7598     0.4444   1.710   0.0873 .  
## Percent_Maple:Percent_Oak   0.5046     1.5930   0.317   0.7514    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Prcn_M Prcn_O
## Percent_Mpl -0.672              
## Percent_Oak -0.614  0.600       
## Prcnt_M:P_O  0.454 -0.691 -0.807
performance::check_overdispersion(model_oak_maple_2)
## # Overdispersion test
## 
##        dispersion ratio =  0.133
##   Pearson's Chi-Squared = 23.914
##                 p-value =      1
#model for Oak, Maple, and Pine together
model_all <- glmer(round(clean_complete) ~ Percent_Maple + Percent_Oak + 
                           Percent_Pine + (1 | trap_name) + (1 | patch_name), 
                         data = stand_ID_filtered, family = poisson)

summary(model_all)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(clean_complete) ~ Percent_Maple + Percent_Oak + Percent_Pine +  
##     (1 | trap_name) + (1 | patch_name)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1860.8    1880.2    -924.4    1848.8       180 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.07728 -0.14342  0.00230  0.09453  0.37783 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  trap_name  (Intercept) 0.5308   0.7286  
##  patch_name (Intercept) 0.2614   0.5112  
## Number of obs: 186, groups:  trap_name, 186; patch_name, 15
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     2.9191     0.3126   9.338  < 2e-16 ***
## Percent_Maple   0.4709     0.5687   0.828   0.4076    
## Percent_Oak     1.4067     0.3558   3.954 7.69e-05 ***
## Percent_Pine    0.9371     0.4249   2.206   0.0274 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Prcn_M Prcn_O
## Percent_Mpl -0.765              
## Percent_Oak -0.741  0.491       
## Percent_Pin -0.744  0.639  0.682
performance::check_overdispersion(model_all)
## # Overdispersion test
## 
##        dispersion ratio =  0.131
##   Pearson's Chi-Squared = 23.652
##                 p-value =      1
## All Variables -----------------------------------------------------------


##Fitting a poisson model with all of the possible response variables, to 
#explore a correlation between all possible variables and moth counts

all_variables_poisson <- glmer(round(clean_complete) ~ (1|trap_name) + 
              Percent_Oak + Percent_Pine + landscape_type + longitude + 
              forest_area_km2 + stand_area_ha + (1|patch_name),
              family = poisson(), data = stand_ID_filtered)

summary(all_variables_poisson)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## round(clean_complete) ~ (1 | trap_name) + Percent_Oak + Percent_Pine +  
##     landscape_type + longitude + forest_area_km2 + stand_area_ha +  
##     (1 | patch_name)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1909.0    1941.5    -944.5    1889.0       180 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.87569 -0.13572  0.00032  0.09462  0.35788 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  trap_name  (Intercept) 0.5393   0.7344  
##  patch_name (Intercept) 0.3123   0.5588  
## Number of obs: 190, groups:  trap_name, 190; patch_name, 15
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)   
## (Intercept)            6.6680214 14.0675971   0.474  0.63550   
## Percent_Oak            0.9611935  0.3182652   3.020  0.00253 **
## Percent_Pine           0.5048426  0.3447402   1.464  0.14308   
## landscape_typeforest   0.9008992  0.3873226   2.326  0.02002 * 
## landscape_typeurban   -0.2759427  0.3971989  -0.695  0.48723   
## longitude              0.0466809  0.1914203   0.244  0.80733   
## forest_area_km2       -0.0005727  0.0003957  -1.447  0.14777   
## stand_area_ha         -0.0134213  0.0193206  -0.695  0.48727   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Prcn_O Prcn_P lndscp_typf lndscp_typr longtd frs__2
## Percent_Oak  0.181                                                    
## Percent_Pin  0.159  0.597                                             
## lndscp_typf -0.236 -0.168 -0.220                                      
## lndscp_typr  0.033  0.139  0.182  0.260                               
## longitude    1.000  0.190  0.166 -0.230       0.044                   
## forst_r_km2  0.493  0.123  0.134 -0.499       0.070       0.497       
## stand_are_h -0.041 -0.048 -0.018 -0.028       0.024      -0.030  0.075
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.243731 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
AIC(all_variables_poisson)
## [1] 1908.988
print(summary(all_variables_poisson))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## round(clean_complete) ~ (1 | trap_name) + Percent_Oak + Percent_Pine +  
##     landscape_type + longitude + forest_area_km2 + stand_area_ha +  
##     (1 | patch_name)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1909.0    1941.5    -944.5    1889.0       180 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.87569 -0.13572  0.00032  0.09462  0.35788 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  trap_name  (Intercept) 0.5393   0.7344  
##  patch_name (Intercept) 0.3123   0.5588  
## Number of obs: 190, groups:  trap_name, 190; patch_name, 15
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)   
## (Intercept)            6.6680214 14.0675971   0.474  0.63550   
## Percent_Oak            0.9611935  0.3182652   3.020  0.00253 **
## Percent_Pine           0.5048426  0.3447402   1.464  0.14308   
## landscape_typeforest   0.9008992  0.3873226   2.326  0.02002 * 
## landscape_typeurban   -0.2759427  0.3971989  -0.695  0.48723   
## longitude              0.0466809  0.1914203   0.244  0.80733   
## forest_area_km2       -0.0005727  0.0003957  -1.447  0.14777   
## stand_area_ha         -0.0134213  0.0193206  -0.695  0.48727   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Prcn_O Prcn_P lndscp_typf lndscp_typr longtd frs__2
## Percent_Oak  0.181                                                    
## Percent_Pin  0.159  0.597                                             
## lndscp_typf -0.236 -0.168 -0.220                                      
## lndscp_typr  0.033  0.139  0.182  0.260                               
## longitude    1.000  0.190  0.166 -0.230       0.044                   
## forst_r_km2  0.493  0.123  0.134 -0.499       0.070       0.497       
## stand_are_h -0.041 -0.048 -0.018 -0.028       0.024      -0.030  0.075
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.243731 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
tab_model(all_variables_poisson,
          show.ci = FALSE,     # hide CI since you only want SE & p
          show.stat = TRUE,    # show test statistics
          p.style = "numeric") # show numeric p-values
  round(clean complete)
Predictors Incidence Rate Ratios Statistic p
(Intercept) 786.84 0.47 0.636
Percent Oak 2.61 3.02 0.003
Percent Pine 1.66 1.46 0.143
landscape typeforest 2.46 2.33 0.020
landscape type [urban] 0.76 -0.69 0.487
longitude 1.05 0.24 0.807
forest area km2 1.00 -1.45 0.148
stand area ha 0.99 -0.69 0.487
Random Effects
σ2 0.56
τ00 trap_name 0.54
τ00 patch_name 0.31
ICC 0.36
N trap_name 190
N patch_name 15
Observations 190
Marginal R2 / Conditional R2 0.211 / 0.495
#checking for collinearity between variables
#using the 'performance' package, which is built for linear 
#mixed models (glmer)
collinearity <- check_collinearity(all_variables_poisson)
print(collinearity)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##             Term  VIF     VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##      Percent_Oak 1.58 [1.35,   1.96]     1.26      0.63     [0.51, 0.74]
##     Percent_Pine 1.65 [1.40,   2.04]     1.28      0.61     [0.49, 0.71]
##   landscape_type 1.57 [1.34,   1.94]     1.25      0.64     [0.51, 0.74]
##        longitude 1.37 [1.20,   1.70]     1.17      0.73     [0.59, 0.84]
##  forest_area_km2 1.78 [1.50,   2.22]     1.33      0.56     [0.45, 0.67]
##    stand_area_ha 1.01 [1.00, 180.64]     1.01      0.99     [0.01, 1.00]
# Convert to a regular data frame
col_df <- as.data.frame(collinearity)
print(col_df)
##              Term      VIF VIF_CI_low VIF_CI_high SE_factor Tolerance
## 1     Percent_Oak 1.581963   1.353069    1.959249  1.257761 0.6321259
## 2    Percent_Pine 1.647135   1.401846    2.042150  1.283408 0.6071148
## 3  landscape_type 1.570228   1.344309    1.944384  1.253087 0.6368502
## 4       longitude 1.370329   1.196894    1.696532  1.170610 0.7297520
## 5 forest_area_km2 1.782123   1.503393    2.215186  1.334962 0.5611286
## 6   stand_area_ha 1.014701   1.000001  180.642302  1.007324 0.9855115
##   Tolerance_CI_low Tolerance_CI_high
## 1      0.510399526         0.7390606
## 2      0.489680063         0.7133451
## 3      0.514301733         0.7438766
## 4      0.589437855         0.8354955
## 5      0.451429289         0.6651623
## 6      0.005535802         0.9999988
# checking for co-variation between numeric predictors
#in a Pearson correlation matrix
numeric_vars <- stand_ID_filtered[, c("Percent_Oak", "Percent_Pine", "longitude", 
                                      "forest_area_km2", "stand_area_ha")]
# Correlation matrix
cor(numeric_vars, use = "complete.obs")
##                 Percent_Oak Percent_Pine   longitude forest_area_km2
## Percent_Oak      1.00000000  -0.52500780 -0.15613141     -0.01062395
## Percent_Pine    -0.52500780   1.00000000 -0.11765413      0.05010221
## longitude       -0.15613141  -0.11765413  1.00000000     -0.41050141
## forest_area_km2 -0.01062395   0.05010221 -0.41050141      1.00000000
## stand_area_ha   -0.03911951   0.09618211  0.03658612     -0.24775945
##                 stand_area_ha
## Percent_Oak       -0.03911951
## Percent_Pine       0.09618211
## longitude          0.03658612
## forest_area_km2   -0.24775945
## stand_area_ha      1.00000000
library(gt)

cor_mat <- cor(numeric_vars, use = "complete.obs")
cor_df <- as.data.frame(round(cor_mat, 3))
cor_df |> gt::gt() |> gt::tab_caption("Correlation Matrix")
Correlation Matrix
Percent_Oak Percent_Pine longitude forest_area_km2 stand_area_ha
1.000 -0.525 -0.156 -0.011 -0.039
-0.525 1.000 -0.118 0.050 0.096
-0.156 -0.118 1.000 -0.411 0.037
-0.011 0.050 -0.411 1.000 -0.248
-0.039 0.096 0.037 -0.248 1.000
#write.csv(cor_df, "correlation_matrix.csv", row.names = TRUE)


#checking for interaction between oak and landscape type
interaction_model <- lm(clean_complete ~ Percent_Oak * landscape_type + 
                          Percent_Pine + longitude + forest_area_km2 + stand_area_ha, 
                        data = stand_ID_filtered)
summary(interaction_model)
## 
## Call:
## lm(formula = clean_complete ~ Percent_Oak * landscape_type + 
##     Percent_Pine + longitude + forest_area_km2 + stand_area_ha, 
##     data = stand_ID_filtered)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -68.68 -32.55 -12.11  22.52 287.54 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                       281.52445  452.00476   0.623  0.53418   
## Percent_Oak                        43.51391   26.86459   1.620  0.10704   
## landscape_typeforest                0.91995   17.39731   0.053  0.95789   
## landscape_typeurban                15.38240   23.83683   0.645  0.51954   
## Percent_Pine                       56.61962   21.26213   2.663  0.00845 **
## longitude                           3.42628    6.17707   0.555  0.57980   
## forest_area_km2                    -0.02668    0.01540  -1.732  0.08493 . 
## stand_area_ha                       0.56936    1.12480   0.506  0.61334   
## Percent_Oak:landscape_typeforest   65.12563   35.49584   1.835  0.06819 . 
## Percent_Oak:landscape_typeurban   -66.97300   54.74304  -1.223  0.22278   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 51.2 on 180 degrees of freedom
##   (55 observations deleted due to missingness)
## Multiple R-squared:  0.1451, Adjusted R-squared:  0.1024 
## F-statistic: 3.396 on 9 and 180 DF,  p-value: 0.0007096
AIC(interaction_model)
## [1] 2046.479
#checking for interaction between oak and landscape type, with a glm model
interaction_model_glm <- glmer(round(clean_complete) ~ 
                                 Percent_Oak * landscape_type + 
                                 Percent_Pine + 
                      (1 | trap_name) + (1 | patch_name), 
                    data = stand_ID_filtered, family = poisson)

summary(interaction_model_glm)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(clean_complete) ~ Percent_Oak * landscape_type + Percent_Pine +  
##     (1 | trap_name) + (1 | patch_name)
##    Data: stand_ID_filtered
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    1898.4    1927.6    -940.2    1880.4       181 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.00781 -0.15339  0.00673  0.09710  0.37024 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  trap_name  (Intercept) 0.5103   0.7143  
##  patch_name (Intercept) 0.3576   0.5980  
## Number of obs: 190, groups:  trap_name, 190; patch_name, 15
## 
## Fixed effects:
##                                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                         2.9314     0.2906  10.089  < 2e-16 ***
## Percent_Oak                         1.1848     0.4106   2.885  0.00391 ** 
## landscape_typeforest                0.5073     0.4099   1.238  0.21587    
## landscape_typeurban                 0.6694     0.5066   1.321  0.18642    
## Percent_Pine                        0.8037     0.3480   2.310  0.02090 *  
## Percent_Oak:landscape_typeforest    0.4032     0.5333   0.756  0.44969    
## Percent_Oak:landscape_typeurban    -2.8278     0.9080  -3.114  0.00184 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                      (Intr) Prcn_O lndscp_typf lndscp_typr Prcn_P
## Percent_Oak          -0.535                                      
## lndscp_typf          -0.526  0.199                               
## lndscp_typr          -0.556  0.328  0.279                        
## Percent_Pin          -0.296  0.365 -0.263       0.230            
## Prcnt_Ok:lndscp_typf  0.268 -0.594 -0.493      -0.125       0.201
## Prcnt_Ok:lndscp_typr  0.244 -0.447 -0.134      -0.573      -0.150
##                      Prcnt_Ok:lndscp_typf
## Percent_Oak                              
## lndscp_typf                              
## lndscp_typr                              
## Percent_Pin                              
## Prcnt_Ok:lndscp_typf                     
## Prcnt_Ok:lndscp_typr  0.267
AIC(interaction_model_glm)
## [1] 1898.357
performance::check_overdispersion(interaction_model_glm)
## # Overdispersion test
## 
##        dispersion ratio =  0.141
##   Pearson's Chi-Squared = 25.499
##                 p-value =      1
performance::check_model(interaction_model_glm)

tab_model(interaction_model_glm,
          show.ci = FALSE,     # hide CI since you only want SE & p
          show.stat = TRUE,    # show test statistics
          p.style = "numeric") # show numeric p-values
  round(clean complete)
Predictors Incidence Rate Ratios Statistic p
(Intercept) 18.75 10.09 <0.001
Percent Oak 3.27 2.89 0.004
landscape typeforest 1.66 1.24 0.216
landscape type [urban] 1.95 1.32 0.186
Percent Pine 2.23 2.31 0.021
Percent Oak × landscape
typeforest
1.50 0.76 0.450
Percent Oak × landscape
type [urban]
0.06 -3.11 0.002
Random Effects
σ2 0.53
τ00 trap_name 0.51
τ00 patch_name 0.36
ICC 0.40
N trap_name 190
N patch_name 15
Observations 190
Marginal R2 / Conditional R2 0.204 / 0.526
ggplot(stand_ID_filtered, aes(x = Percent_Oak, y = clean_complete, color = landscape_type)) +
  geom_point(alpha = 0.6) +  # Add raw data points with transparency
  geom_smooth(method = "lm", se = FALSE) +  # Add regression lines by group
  theme_minimal() +
  labs(
    title = " ",
    x = "Percent Oak",
    y = "Moth Counts",
    color = "Landscape Type"
  ) +
  scale_color_viridis_d(option = "D")  # You can also try options "A", "B", "C", "E"

# Heatmap oak/pine spectrum -----------------------------------------------

# Create the 'oak/pine spectrum' table
table(stand_category_filtered$Percent_Oak, 
      stand_category_filtered$Percent_Pine)
##      
##        0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
##   0    0   0   0   0   1   4   4   8   6
##   0.1  0   0   0   0   0   2   3   7   0
##   0.2  0   0   0   1   5   1   0   0   0
##   0.3  0   0  10   6  12  12   2   0   0
##   0.4  5  12  12   5   1   0   0   0   0
##   0.5 14   2   0   0   2   0   0   0   0
##   0.6 15   3   0   0   0   0   0   0   0
##   0.7 12   7   0   0   0   0   0   0   0
##   0.8 10   2   0   0   0   0   0   0   0
# Create the contingency table
contingency_table <- table(stand_type_filtered$Percent_Oak, 
                           stand_type_filtered$Percent_Pine)

# Convert the table to a data frame
contingency_df <- as.data.frame(contingency_table)

# Create a plot (heatmap) of the contingency table
contingency_df 
##    Var1 Var2 Freq
## 1     0    0   19
## 2   0.1    0    4
## 3   0.2    0    2
## 4   0.3    0   13
## 5   0.4    0   20
## 6   0.5    0   14
## 7   0.6    0   15
## 8   0.7    0   12
## 9   0.8    0   10
## 10    0  0.1    0
## 11  0.1  0.1    2
## 12  0.2  0.1    0
## 13  0.3  0.1    0
## 14  0.4  0.1   12
## 15  0.5  0.1    2
## 16  0.6  0.1    3
## 17  0.7  0.1    7
## 18  0.8  0.1    2
## 19    0  0.2    3
## 20  0.1  0.2    0
## 21  0.2  0.2    0
## 22  0.3  0.2   11
## 23  0.4  0.2   12
## 24  0.5  0.2    0
## 25  0.6  0.2    0
## 26  0.7  0.2    0
## 27  0.8  0.2    0
## 28    0  0.3    0
## 29  0.1  0.3    0
## 30  0.2  0.3    1
## 31  0.3  0.3    6
## 32  0.4  0.3    5
## 33  0.5  0.3    0
## 34  0.6  0.3    0
## 35  0.7  0.3    0
## 36  0.8  0.3    0
## 37    0  0.4    1
## 38  0.1  0.4    0
## 39  0.2  0.4    5
## 40  0.3  0.4   12
## 41  0.4  0.4    1
## 42  0.5  0.4    2
## 43  0.6  0.4    0
## 44  0.7  0.4    0
## 45  0.8  0.4    0
## 46    0  0.5    4
## 47  0.1  0.5    2
## 48  0.2  0.5    1
## 49  0.3  0.5   12
## 50  0.4  0.5    0
## 51  0.5  0.5    0
## 52  0.6  0.5    0
## 53  0.7  0.5    0
## 54  0.8  0.5    0
## 55    0  0.6    4
## 56  0.1  0.6    3
## 57  0.2  0.6    0
## 58  0.3  0.6    2
## 59  0.4  0.6    0
## 60  0.5  0.6    0
## 61  0.6  0.6    0
## 62  0.7  0.6    0
## 63  0.8  0.6    0
## 64    0  0.7    8
## 65  0.1  0.7    7
## 66  0.2  0.7    0
## 67  0.3  0.7    0
## 68  0.4  0.7    0
## 69  0.5  0.7    0
## 70  0.6  0.7    0
## 71  0.7  0.7    0
## 72  0.8  0.7    0
## 73    0  0.8    6
## 74  0.1  0.8    0
## 75  0.2  0.8    0
## 76  0.3  0.8    0
## 77  0.4  0.8    0
## 78  0.5  0.8    0
## 79  0.6  0.8    0
## 80  0.7  0.8    0
## 81  0.8  0.8    0
ggplot(contingency_df, aes(x = Var1, y = Var2, fill = Freq)) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "darkred") +
  labs(x = "Proportion Oak", y = "Proportion Pine", fill = "Frequency") +
  theme_minimal()

# Save the plot as an image (e.g., PNG)
#ggsave("oak_pine_heatmap.png", width = 7.5, height = 6)


# Citations ---------------------------------------------------------------

##   @Manual{,
##     title = {R: A Language and Environment for Statistical Computing},
##     author = {{R Core Team}},
##     organization = {R Foundation for Statistical Computing},
##     address = {Vienna, Austria},
##     year = {2021},
##     url = {https://www.R-project.org/},
##   }

citation("lme4")
## To cite lme4 in publications use:
## 
##   Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015).
##   Fitting Linear Mixed-Effects Models Using lme4. Journal of
##   Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Article{,
##     title = {Fitting Linear Mixed-Effects Models Using {lme4}},
##     author = {Douglas Bates and Martin M{\"a}chler and Ben Bolker and Steve Walker},
##     journal = {Journal of Statistical Software},
##     year = {2015},
##     volume = {67},
##     number = {1},
##     pages = {1--48},
##     doi = {10.18637/jss.v067.i01},
##   }

xfun::Rscript_call( rmarkdown::render, list( input = “Report/2023_2024_All_Data_Oct_25.Rmd”, output_format = “html_document” ), fail = TRUE )

xfun::Rscript_call( rmarkdown::render, list(input = “Report/2023_2024_All_Data_Oct_25.Rmd”, output_format = “html_document”) )

xfun::Rscript_call( rmarkdown::render, list(input = “Report/2023_2024_All_Data_Oct_25.Rmd”, output_format = “pdf_document”) )